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We study critical random Boolean networks with two inputs per node that contain only canalyzing 
functions. We present a phenomenological theory that explains how a frozen core of nodes that are 
frozen on all attractors arises. This theory leads to an intuitive understanding of the system's 
dynamics as it demonstrates the analogy between standard random Boolean networks and networks 
with canalyzing functions only. It reproduces correctly the scaling of the number of nonfrozen nodes 
with system size. We then investigate numerically attractor lengths and numbers, and explain the 
findings in terms of the properties of relevant components. In particular we show that canalyzing 
networks can contain very long attractors, albeit they occur less often than in standard networks. 
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1. INTRODUCTION 

Random Boolean networks are often used as generic 
models for the dynamics of complex systems of interact- 
ing entities, such as social and economic networks, neural 
networks, and gene or protein interaction networks [1]. 
The simplest and most widely studied of these models 
was introduced in 1969 by Kauffman [2] as a model for 
gene regulation. The system consists of N nodes, each 
of which receives input from K randomly chosen other 
nodes. The network is updated synchronously, the state 
of a node at time step t being a Boolean function of the 
states of the K input nodes at the previous time step, 
t — 1. The Boolean update functions are randomly as- 
signed to every node in the network, and together with 
the connectivity pattern they define the realization of the 
network. For any initial condition, the network eventu- 
ally settles on a periodic attractor. Of special interest 
are critical networks, which lie at the boundary between 
a frozen phase and a chaotic phase [3-5]. In the frozen 
phase, a perturbation at one node propagates during one 
time step on an average to less than one node, and the at- 
tractor lengths remain finite in the limit N — > oo. In the 
chaotic phase, the difference between two almost identi- 
cal states increases exponentially fast, because a pertur- 
bation propagates on an average to more than one node 
during one time step [6] . 

Critical networks with K = 2 inputs per node have 
been studied by a variety of authors. A K = 2 net- 
work is critical if frozen and reversible update functions 
are chosen with equal probability. The remaining up- 
date functions are canalyzing, i.e., one input can fix the 
output of a node, irrespective of the value of the second 
input. Table I shows the 16 update functions of K = 2 
networks. Critical networks that contain a nonvanishing 
proportion of frozen and reversible update functions are 
in the meantime relatively well understood, see [7-13]. 
They contain three classes of nodes, which behave differ- 
ently on attractors. First, there are nodes that are frozen 
on the same value on every attractor. Such nodes give 
a constant input to other nodes and are otherwise irrel- 
evant. They form the frozen core of the network. Sec- 
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TABLE I: The 16 update functions for nodes with two in- 
puts. The first column lists the 4 possible states of the two 
inputs, the other columns represent one update function each, 
falling into the classes frozen {J-), canalyzing (Ci and C2) and 
reversible (1Z)- 



ond, there are nodes without outputs and nodes whose 
outputs go only to irrelevant nodes. Though they may 
fluctuate, they are also classified as irrelevant since they 
act only as slaves to the nodes determining the attractor 
period. Third, the relevant nodes are the nodes whose 
state is not constant and that control at least one rele- 
vant node. These nodes determine completely the num- 
ber and period of attractors. If only these nodes and 
the links between them are considered, these nodes form 
loops with possibly additional links and chains of rele- 
vant nodes within and between loops. We call a set of 
relevant nodes that are connected in this way a relevant 
component. The nonfrozen nodes that are not relevant 
sit on trees rooted in the relevant components. In [8], it 
was found that the number of nonfrozen nodes scales in 
the limit N — ► 00 as TV 2 / 3 and the number of relevant 
nodes as N 1 / 3 . This result was confirmed by an analyt- 
ical calculation in [13], where it was additionally shown 
that the number of nonfrozen nodes with two nonfrozen 
inputs scales as N 1 ^ 3 , and that the number of relevant 
nodes with two relevant inputs remains finite in the limit 
A" — > 00. The mean number of relevant components was 
found to be proportional to In N, and all but the largest 
relevant components are simple loops. 

Canalyzing networks share many features of other crit- 
ical networks. Thus, the calculation by Samuelsson and 
Troein [9] of the number of attractors can be general- 
ized to canalyzing networks [12], implying that canalyz- 
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ing networks also have of the order of TV 2 / 3 nonfrozcn 
nodes and at most TV 1 / 3 nonfrozen nodes with two rel- 
evant inputs, and that the mean attractor number in- 
creases faster than any power law with the network size. 
Whether the nonfrozen nodes are the same on all attrac- 
tors in canalyzing networks (as is the case of networks 
containing frozen update functions), can not be decided 
from previous work. In particular, the detailed results of 
[13] can only be derived if there are nodes with frozen 
functions. For these reasons, a separate study of can- 
alyzing networks is needed. It is the main aim of this 
paper to show how the attractors and the frozen nodes 
arise in canalyzing networks. We will see that canalyzing 
networks also have a frozen core, which means that most 
frozen nodes are the same on all attractors. It follows 
then that all results of [13] about the relevant part of the 
network can be applied also to canalyzing networks. We 
will also put an end to the long-standing belief that can- 
alyzing networks have less and shorter attractors. These 
features were argued to make canalyzing networks bio- 
logically more relevant [1]. 

Let us therefore focus on K = 2 networks that contain 
only C2 functions. These functions take one value (0 or 
1) three times and the other one once. This means that 
each of the two inputs can fix the output of the function 
irrespective of the other input. For instance, the out- 
put of the first C2 function shown in Table I must be 
if the first input is 1, irrespective of the second input. 
It must also be if the second input is 1, irrespective 
of the first input. Each of the 8 C2 functions is chosen 
with equal probability in our simulations. We will com- 
pare our results for C2 networks with those of standard 
random Boolean networks (RBNs), where all 16 update 
functions have the same weight. Part of our results will 
also be compared to those of C\ networks, where the up- 
date functions are chosen only from the C\ class. The 
Ci networks can be trivially mapped on critical K = 1 
networks by removing the link to the input to which the 
node does not respond. These networks have no frozen 
core. They have of the order of VN relevant nodes, ar- 
ranged in ~ In N simple loops (see [14]), with the largest 
loop length being of the order >/~N. The other nodes sit 
on trees rooted in these loops. 

In the next Section, we study numerically the frozen 
nodes in order to find out if the same nodes are frozen on 
all attractors of C2 networks. In Section 3, we explain the 
results of the numerical simulations using phenomenolog- 
ical arguments and analytical calculations. In Section 4, 
we study the number and length of attractors of C2 net- 
works and compare the results to those of other network 
types. Finally we summarize and discuss our results in 
the last Section. 



2. THE FROZEN CORE 

From a generalization of the work of Samuelsson and 
Troein [9] to all critical K = 2 random Boolean networks 




FIG. 1: Mean number of nonfrozen nodes for canalyzing C2 
networks (stars) and RBNs (circles). Solid lines connect the 
data points for N—Nf (mean number of nonfrozen nodes per 

attractor), dashed lines the data points for N — iVl (mean 
number of nodes that are nonfrozen at least on one found 
attractor). The dotted line is a power law N 2/:i . Different 
attractors are counted only once, without considering their 
basins of attraction. We have considered 1000 initial states 
per network, averaged over 2000 networks. 

[12], we know that for canalyzing networks the number of 
nonfrozen nodes scales for large network size in the same 
way as for RBNs, i.e., with TV 2 / 3 . In RBNs, the nodes 
frozen on all attractors (i.e., the nodes belonging to the 
frozen core) can be identified by starting with the nodes 
with frozen update functions and by iteratively determin- 
ing nodes that become frozen because of frozen inputs, 
see [13]. In canalyzing networks there are no frozen func- 
tions to start with, and this method cannot be applied. It 
therefore arises the question whether canalyzing networks 
have a frozen core at all, or whether different attractors 
have different nonfrozen nodes. In the following, we will 
answer this question using computer simulations. 

In Fig. 1 we show the average number of frozen 
nodes per attractor, both for canalyzing networks and 
for RBNs. We actually plot the difference (N - N ( f a) ) 
as a function of N in order to better see the asymptotic 
behavior. The other two curves show the difference (N — 
Nj), where Nj is the number of nodes frozen on all 
attractors found in the simulation of a network. The 
technical details can be found in the caption to the figure. 

One can clearly see the similarity of the results for the 
mean number of nonfrozen nodes per attractor (N—Nf) 
for canalyzing networks and for RBN, in agreement with 
[9, 12]. The expected power law TV 2 / 3 is not yet reached 
for the system sizes used in the simulation and is only 
approached slowly with increasing N. The number of 
nonfrozen nodes in the simulations was only of the order 
of 100 for the largest simulated networks, which is yet too 
small to see the asymptotic behavior. (And the number 
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of relevant nodes, which increases as JV 1 / 3 , is only of the 
order 10 for the largest simulated networks.) 

Our results for (JV — JVi ) suggest that canalyzing net- 
works have a frozen core and of the order of JV 2 / 3 non- 
frozen nodes, because the curves for (JV - N { f n} ) differ 

only by a constant factor from those for (JV — JVj ) for 
both network types. Furthermore, there are of the or- 
der of TV 2 / 3 nodes that are frozen only on part of the 
attractors. The factor between the two curves is larger 
for canalyzing networks than for RBN. 

In the following, we will explain the reason for the 
constant factor between the curves for (JV - N { f a} ) and 

(JV - JVj n) ). Since this point has not yet been discussed 
for RBNs, we will consider both RBNs and canalyzing 
networks. The explanation of the origin of the frozen 
core in canalyzing networks will be postponed until the 
next Section. 

The difference between the curves for (JV- Jvj a) ) and 

(JV — Nj) is due to those nodes that are frozen on some 
attractors, but not on all attractors. These nodes do not 
belong to the frozen core, and they are therefore rele- 
vant nodes or sit on nonfrozen trees that are rooted in 
relevant components. In Section 1, we have mentioned 
that most relevant components consist of simple loops, 
and that only a few large components arc more complex 
and contain relevant nodes that have two relevant in- 
puts. Clearly, since the dynamics of the nonfrozen nodes 
in the trees is determined by the dynamics of the rele- 
vant nodes, all nodes of a relevant component and the 
nonfrozen trees rooted in it undergo a cycle of the same 
period (when the attractor has been reached), which is 
determined by the initial state of the relevant nodes of 
that component. If this cycle has period one, all nodes of 
this component are frozen on this attractor. We there- 
fore have to show that a finite fraction of all nonfrozen 
nodes are on cycles of length 1 on a finite fraction of all 
attractors. This would lead to a constant factor between 
the curves for (JV - N { f a} ) and (JV - Jvj n) ). 

Let us first consider relevant components that are sim- 
ple loops, and their nonfrozen trees. The mean num- 
ber of relevant loops of length I is 1/7 for all I up to a 
cutoff l c ~ JV 1 / 3 . The mean size of a tree rooted in a 
relevant node is JV 1 / 3 . The largest of these components 
consists therefore of the order of JV 2 / 3 nodes (including 
the nonfrozen trees), and if such a component reaches a 
fixed point attractor with nonzero probability, we have 
explained the factor between the two curves. A relevant 
loop of length I has either two fixed points (if the loop 
is "even", i.e. if the state of a node is repeated after I 
time steps) or none (if the loop is "odd" , i.e. the state 
of a node is inverted after I time steps). Each case oc- 
curs with probability 1/2. The number of attractors of 
a component with a loop of length I, however, increases 
exponentially with /, and for this reason only a vanishing 
proportion of attractors of components of a size of the 
order l c are fixed points in the limit of large JV. 
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FIG. 2: Number of frozen nodes for K = 1 networks. Large 
squares: iVy , small squares: JVy . For comparison, the cor- 
responding curves for RBNs and C2 networks are also shown, 
the data is taken from Fig. 1. The dotted lines correspond to 
the power laws JV, \/N and \/N /2. 



Next, we consider complex relevant components. In 
contrast to simple loops, where each initial state is part 
of a periodic cycle in state space, more complex compo- 
nents can have fixed points that are true attractors, i.e. 
that are reached from a nonvanishing proportion of ini- 
tial states (but not from all initial states). One example 
of such a component in RBNs was discussed in [15]. It 
is a loop with an additional chain of nodes within the 
loop, such that there is one node that has two relevant 
inputs. From [13], we know that this component occurs 
with nonvanishing probability in an RBN. In the case 
that the update function of the node with two inputs is 
only if both inputs are and that the two numbers of 
nodes on the two subloops have a common divisor greater 
1, all apart from a finite number of initial conditions end 
up on the same fixed point. The existence of such com- 
ponents docs not only explain the multiplicative factor 
between the curves for (JV - Jvj n) ) and (JV - Jvj a) ), but 
it explains also why this factor is larger for C2 networks 
than for RBNs. The probability that an update function 
as in the above example is chosen at the relevant node 
with two inputs is larger for canalyzing networks. 

We conclude this Section by comparing C2 networks 
with C\ networks, which are canalyzing networks, but 
with update functions of the C\ class. Our simulation re- 
sults for JVj Q) and Jvj n) are shown in Fig. 2. Both curves 

increase for C\ networks as s/N. Asymptotically, JV^ is 

of the order VJV /2, since only even loops of length 1 can 
be frozen on all attractors, while the average tree size 
is of the order of \/JV. To calculate Jvj°' ) we note that 
only even loops of length I have (two) fixed points, and 
that they do not reach these fixed points for 2 l — 2 initial 
conditions. The contribution to the average number of 
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frozen nodes per attractor from simple loops of length I 
is then l/l ■ 1/2 • 21 /2 l ■ s/N, where we again take the 
probability of an even loop 1/2 and the average size of 
a tree of the order of s/N into account. Summing over I 
leads to Nj — > t/N. Wc conclude that, differently from 

C 2 networks, the ratio (N - N { f n) )/(N - M o) ) approaches 
1 asymptotically. As we have learned, the reason for a 
larger factor between the two curves for RBNs and C2 net- 
works is the existence of complex relevant components. 
In C\ networks all relevant components are simple loops. 



3. SELF-FREEZING LOOPS 

In this Section, we want to explain how the frozen core 
arises in Ci networks and find some of its properties. We 
also estimate its size by means of analytical arguments. 
The results are in agreement with [9, 12], confirming our 
intuitive understanding of the origin of the frozen core. 

Since a C2 network has no nodes with a frozen function, 
the frozen core cannot be formed starting from single 
frozen nodes. Instead, there must exist groups of nodes 
that fix each other's value and do not respond to changes 
in nodes outside this group. Let us consider the simplest 
example of such a group, namely a loop where each node 
canalyzes (fixes) the state of its successor once it settles 
on its majority bit (the one occuring 3 times in its update 
function table). We call such loops self- freezing loops. In 
the following, we first discuss these self-freezing loops, 
before discussing how a frozen core that contains almost 
all nodes can be formed starting from these loops. 

If all nodes of a self-freezing loop are on their majority 
bits, it stays frozen forever. Starting from an arbitrary 
initial state, the number of nodes of a self-freezing loop on 
majority bits cannot decrease with time, since each such 
node drives its successor to its majority bit. It remains 
constant only in the unlikely case that all inputs from 
outside the loop to the nodes of the loop are fixed on 
the noncanalyzing value. We can therefore assume that 
self- freezing loops are usually frozen on all attractors, at 
least if the loops are large. As we will see, most nodes 
that are part of self-freezing loops sit in loops with a size 
of the order of A 1 / 3 . 

The number of nodes on self-freezing loops can be es- 
timated as follows. The probability that a given node 
constitutes a self-freezing loop of length 1 is 1/N for a 
network with A nodes. It is the product of the probabil- 
ity 2/N that the node is self-connected and the proba- 
bility 1/2 that the node is canalyzed by its own majority 
bit. There is thus on average one self- freezing loop of 
length 1 per network. With the same line of reasoning, 
the average number of self-freezing loops of length 2 per 
network is obtained to be ( 2 ) (5) ~ \- For the 

self-freezing loops of length I > 2 one has to take into 
account an additional factor, corresponding to the num- 
ber of ways to construct a directed loop out of I nodes. 
The number of self-freezing loops of length I per network 



is found to be l/l. The overall number of nodes on self- 
freezing loops /o is then 

f0=4=l\l = lc (1) 

l c being the cutoff in loop length. This simple proba- 
bilistic considerations applies if l c is much smaller than 
N. 

We can obtain a confirmation of this estimation and a 
result for the value of l c by mapping the problem of find- 
ing a self-freezing loop in a C2 network onto the problem 
of finding the relevant nodes sitting on relevant loops in 
a critical network that contains no canalyzing functions 
at all, but only 1Z and T functions. Whether a randomly 
chosen node in such a network is part of a relevant loop 
is determined by the following algorithm. Consider the 
two inputs of this node. With probability 1/4, both in- 
puts have a frozen update function, and the node is not 
relevant. With probability 1/2, one input has a frozen 
update function and the other one a reversible one. In 
this case we draw a link to this reversible input node 
and thus mark it for investigation of its two inputs in 
the next step. With probability 1/4, both inputs have 
reversible update functions, and we draw links to both of 
them. We iterate this procedure, choosing at each step 
the two inputs of a node at random from all nodes, and 
drawing links to those inputs that do not have frozen up- 
date function. The procedure continues until we either 
find a connection back to the original node (in which 
case it is relevant), or until no more links can be drawn 
(in which case the original node is not part of a relevant 
loop). Prom the results of our article [13], we know that 
there is a mean number of l/l relevant loops of size I in 
such a network, and that the cutoff in the size of relevant 
components scales as A 1 / 3 . 

Now, we turn to the procedure of finding self-freezing 
loops in C2 networks and show that it is identical to 
the procedure just described. We start with a randomly 
chosen node and determine whether it is part of a self- 
freezing loop. Consider the two inputs of this node. With 
probability 1/4, the majority bit of neither input cana- 
lyzes the chosen node, and the node is not on a self- 
freezing loop. With probability 1/2, the majority bit of 
one input does not canalyze the chosen node, and the 
majority bit of the other input canalyzes it. Let us draw 
a link to this input node and consider its two inputs in 
the next step. With probability 1/4, the majority bits of 
both inputs canalyze the chosen node, and we draw links 
to both of them. We iterate this procedure, choosing at 
each step the two inputs of a node at random from all 
nodes, and drawing links to those inputs, whose major- 
ity bits canalyze the node. The procedure continues until 
we either find a connection back to the original node (in 
which case it is part of a self-freezing loop), or until no 
more links can be drawn (in which case the original node 
is not part of a self-freezing loop). The analogy of the 
two procedures is obvious, and we conclude l c ~ A 1 / 3 
and f ~ A 1 / 3 . 
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Obviously, nodes depending on or canalyzed by the 
frozen nodes of the self-freezing loops freeze also, and 
such nodes may lead to the freezing of further nodes, etc. 
We introduce a dynamical process in order to determine 
the total number of nodes that become frozen because of 
the self-freezing loops. We denote with / the number of 
nodes that have already become frozen during the pro- 
cess, and the influence of which on other nodes has yet 
to be determined, ni is the number of nodes for which 
we already know that one of their inputs is frozen but 
does not canalyze them, and 712 is the number of nodes 
for which no frozen input was yet identified. Initially, 
ni = 0, / = fa, n 2 = N - f and n = f + n x + n 2 = N. 
We answer for one of the frozen nodes at a time the 
question whose input it is. It is an input to any of the 
n 2 nodes with probability 2/n. With equal probability 
1/2, the node either becomes frozen by this input, or 
it becomes a non-frozen node with effectively one input 
(called a C\ node in the following). If a node with one 
input chooses the given frozen node as input (that hap- 
pens with probability 1/rt), it becomes frozen. At each 
step, the connected frozen node is being excluded from 
further consideration. The dynamical process stops when 
all the nodes are frozen (which is improbable) or when 
there are no more frozen nodes the influence of which on 
other nodes has not yet been determined. The dynamical 
equations for this process are 

Am = (n 2 - nx)/n, (2) 
/S.1%2 = — 2riiln . 

The stochastic Poisson distributed noise terms £1 and £2 
with the mean values n\jn and n2/n respectively have 
to be taken into account in the equation for /, since / 
becomes small during the process, so that the noise be- 
comes important, compare [13]. The sum n = f+ni+n 2 
decreases by 1 in each step. 

Simulations of this process show that the total number 
of nodes that are frozen because of the self-freezing loops 
is around ~ A 0,8 , and that the number or nodes that are 
not fixed by the self-freezing loops is of the order of N. 
The number of nodes frozen because of the self-freezing 
loops is not large enough to explain the simulation data of 
the previous Section. We therefore have to find a mecha- 
nism that generates more frozen nodes. This is found by 
extending the definition of self-freezing loops. We have 
just seen that nodes with one nonfrozen input appear as 
we identify frozen nodes. Among the nodes that are not 
frozen by the original self-freezing loops, there are new 
types of self-freezing loops that contain chains of nodes 
with one nonfrozen input between C2 nodes. If a chain 
between two C2 nodes as a whole inverts its input, the 
inverted majority bit of the first C2 node has to canalyze 
the second C 2 node. As with original self-freezing loops 
we can claim that the generalized self-freezing loops are 
usually frozen on attractors. At the end of the process 
described by (2), the generalized self- freezing loops need 



to be found. The only effect of nodes with C\ functions 
in such loops is to delay the signal propagation between 
two adjacent C2 nodes. The remaining 712 nodes with 
C2 functions can therefore be considered as an effective 
C2 network, which leads to 712 1 / 3 nodes on generalized 
self-freezing loops with similar loop size statistics as dis- 
cussed above. The independence of the second search for 
self-freezing loops from the first one is due to the large 
enough number m at the end of the process (2). This n\ 
ensures that typical self-freezing loop in the effective C 2 
network have insertions of C\ chains. 

With the new self-freezing loops we again run the dy- 
namical process of type (2). We can even take over the 
values of m, ri2 and n = Hi + n,2 at the end of the first 
process, since tt.2 1 / 3 frozen nodes are negligible in com- 
parison with n,2. Therefore the two equations 

Ani = (ri2 - n{)/n, 

Ari2 = — 2ri2/n. (3) 

apply to both processes together. The equation for n is 
An = — 1, as before. The solution of these equations 
is obtained by going to differential equations for dni/dn 
and dn2/dn, which have the solution 



m = n - — . (5) 

In the same way, at the end of the second process we 
have again an effective C2 network, with chains containing 
newly generated C\ nodes. The number of remaining C\ 
nodes increases in the second process, the number of C2 
nodes decreases, thus leading to an increasing weight of 
C\ nodes in the nonfrozen network. Equations (3) are 
now applied to a third process, similar to the first two. 

The repeated process of identifying generalized self- 
freezing loops and the nodes frozen by them breaks down 
when the remaining nonfrozen nodes cannot be consid- 
ered as an effective C2 network any more. This happens 
when the proportion 712/n of C2 nodes becomes of the 
order ~ 1/y/n. Then, in the process of building a self- 
freezing loop, there occur C\ chains of an average length 
of the order ~ ^fn between C2 nodes. Now, the proba- 
bility to attach a C2 input at the end of the chain is of 
the same order of magnitude 1 / yfn as the probability to 
attach some node of this chain at the end of the chain, in 
which case the chain becomes a loop, and the assembly 
of self-freezing loop becomes improbable. 

Let us denote by N n f the average number of nonfrozen 
nodes in C2 networks and by N2 the average number of 
nonfrozen nodes with two inputs. Ni = N n f — N2 is 
then the average number of nodes with one nonfrozen 
input. The breakdown condition for the iterated process 
becomes then N2 ~ -J 'N n f. Inserting the condition n. 2 ~ 
\J~n in the solution (4) , we obtain 

N n f ~ A 2 / 3 , (6) 
N 2 ~ A 1/3 . 
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This is in agreement with the results of [9, 12] and con- 
firms our intuitive understanding of the frozen core. The 
frozen core consists of sclf-frcczing loops, which arise in 
the interated process described in this Section. The num- 
ber of nonfrozen nodes is of the order TV 2 ' 3 , with only 
TV 1 / 3 nonfrozen nodes having 2 nonfrozen inputs. The 
properties of the nonfrozen part of the network are there- 
fore the same as those of RBNs, and we can take over 
the results obtained for the nonfrozen part of RBNs. 
In particular, we can conclude that the number of rel- 
evant nodes scales as TV 1 / 3 , with only a finite number of 
them having two nonfrozen inputs, and with most rele- 
vant components being simple loops. 

The considerations of this Section can be repeated 
without change also for mixed C\ and C2 critical Boolean 
networks, consisting of nonfrozen nodes with one input 
and of nodes with two inputs having update functions of 
class C2 , provided that the number of nodes with two 
inputs is larger than y/~N (otherwise we are left with a C\ 
network). Therefore, all the results valid for C2 networks 
apply also to mixed C\ /C2 networks. 



4. NUMBER AND LENGTH OF ATTRACTORS 

We also performed simulations to obtain statistical 
properties of the attractors of C2 networks in compari- 
son to RBNs and C\ networks. With the intuitive un- 
derstanding developed in the previous Section we can 
interpret the results and gain some new insights. 

We start with probability density for the attractor 
lengths. The results are presented in Fig. 3. In order to 
understand them we remind ourselves that the number 
of relevant nodes N re i scales in C2 networks and RBNs 
with TV 1 / 3 , whereas for C\ networks it scales with y/N. 
In all cases relevant components of sizes less than ~ iV re | 
are mainly simple loops. The mean number of relevant 
loops of length Z is 1/7 as long as I is sufficiently far below 
N re i . In the last graph of Fig. 3 (lower right corner) the 
curves for TV = 512 for RBNs and C2 networks agree well 
with the curve for C\ networks for N = 64 for smaller 
L. The reason is that the number of relevant nodes is 
for all three types of networks of the same size (since 
512 1 / 3 = 64 1 / 2 ). The difference for large L is due to the 
fact that in C\ networks all relevant nodes are on simple 
loops, so that all relevant components have a cycle period 
of the order of their size, while the more complex compo- 
nents occuring in the other network types can have much 
longer cycle lengths. For large L, the small difference be- 
tween C2 networks and RBNs is due to the presence of 
reversible functions in RBNs: relevant components con- 
taining relevant nodes with two relevant inputs and a 
reversible update function can have extremely long peri- 
ods, which can become of the order of the state space of 
the component [15]. 

In our simulations, we find very long attractors also 
for C2 networks. This fact is fairly surprising if we re- 
member that C2 networks were originally thought to be 
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FIG. 3: Probability density distribution for the attractor 
lengths at different fixed network sizes N. The four figures 
correspond respectively to Ci networks, RBNs, canalyzing C2 
networks and the comparison of the three classes of networks. 
In the lower right figure Ci networks with 64 nodes have on 
average the same number of relevant nodes as C2 networks 
with 512 nodes. In calculating the relative frequencies of at- 
tractor lengths, the frequencies were weighted with the sizes 
of attractors' basins of attraction. We considered 1000 ini- 
tial states for one network realization and averaged over 2000 
networks. The data is binned on a logarithmic scale. 



interesting for their short attractors, see f.i. [2]. We can 
explain the appearance of very long attractors in can- 
alyzing networks in the following way. Let us consider 
a relevant component of an RBN that is a loop with 
an additional chain of nodes within the loop and a re- 
versible update function at the node with two inputs. As 
shown in [15], the attractors of such a component can 
comprise a finite proportion of the state space even for 
very large components. Now, it is possible to construct a 
relevant component that contains three nodes with two 
relevant inputs, which all have a canalyzing function, and 
that has, for the mapping implied by Fig. 4, exactly the 
same attractor states as the relevant component of the 
RBN. The number and lengths of attractors is therefore 
identical in the two components. The reversible function 
constructed from three canalyzing functions does not, of 
course, appear as often in C2 networks as a reversible 
function occurs in the RBN. Therefore the very long at- 
tractors appear relatively seldom in canalyzing networks. 

In order to produce the curves in Fig. 3 we have 
binned the data on a logarithmic scale, and we have cho- 
sen the binning interval large enough to smoothen quite 
large fluctuations. Otherwise, we could see even-odd fluc- 
tuations in the number of attractors with neighboring 
lengths, see Fig. 5. There are always more attractors 
with even lengths. Let us explain this behavior. The 
small components of the relevant part of C2 networks are 
simple loops. But for simple loops even attractor lengths 
appear more often, since a loop with TV; nodes leads to 
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FIG. 4: An RBN relevant component and a C2 relevant com- 
ponent, whose attractors can be mapped pairwise onto each 
other. Triangle arrows represent a chain of nodes, and without 
loss of generality we can assume that the update functions of 
the nodes with one input are all a "copy" function. The node 
marked with its "copy" function is absent in the correspond- 
ing right arch of the component on the right, and that arch 
is thus shorter by one node. The left arch and the straight 
chain of the left component are identical to the two left arches 
and the two straight chains of the right component. The table 
explains the notation for the update functions, 9 being the re- 
versible one, which can be emulated by using the canalyzing 
functions 1,8, and 14. At the node with the update function 
14 the binary input combination 11 never occurs. 




FIG. 5: Small-L part of the graphs for C2 networks (left) and 
RBNs (right) of Fig. 3, without data binning, making visible 
even-odd oscillations. 



FIG. 6: Mean attractor length for canalyzing C2 networks 
(stars), RBNs (circles) and Ci networks (squares) as a func- 
tion of network size N. The dotted line is a power law i/~N . 
For the calculation of the mean attractor length, the attrac- 
tor lengths of different attractors (obtained with 1000 initial 
states for 2000 networks) were weighted with the correspond- 
ing basins of attraction. 



an overestimation of the average attractor length for the 
considered value of N , compare Fig. 6 below. 

Fig. 6 shows our results for the mean length of attrac- 
tors. The most important observation is that the mean 
attractor length grows faster than any power law with 
N for all the considered network classes, including C2 
networks, whose attractors are not substantially shorter 
than those of the other networks. Only for small N can 
one roughly fit the curves with the y/~N law suggested a 
long time ago. 

All the curves are similar in shape. The reason for 
this is again that the relevant part of all three network 
types consists mainly of simple loops. The mean attrac- 
tor length of RBNs becomes larger than that of the other 
two network classes for large N, due to the reversible 
functions occuring in complex components. 

We now discuss simulation results for attractor num- 
bers. In analytical calculations [9, 12], the average num- 
ber CL{N re i) of attractors of length L was found to scale 
with a power of the number of relevant nodes, 



C L (N ra )~N, 



re! 



(7) 



attractors of length N or 2Ni , depending on whether it 
is even or odd. Therefore, an even attractor length 2Ni 
occurs in loops with 2Ni or A"; nodes. 

The cutoff in observed attractor length is of the or- 
der L ~ AN 2 , with A = 0.1 for RBN and 0.01 for C 2 
networks. This is a finite size effect. Also, full ensemble 
averages are hard to reproduce numerically. For example, 
with a substantial probability some network realizations 
appear with untypically large attractors, which lead to 



with a proportionality factor that depends on L. Gl is 
closely related to the number of different possible cycles 
of length L on simple loops. One can approximately write 
Gl ~ 2 L /L, for details see [12], at the end of Section 2. 

Fig. 7 shows our simulation results for the number 
of attractors of length L = 7 and L — 8 for C2 and C\ 
networks. The data do not contradict the theoretically 
predicted scaling with N re i (7), if one realizes the limita- 
tions of computer simulations. We considered 1000 ini- 
tial states and could therefore explore the state space of 
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FIG. 7: Absolute number of attractors of the lengths L = 7 
(dashed and dotted lines) and L = 8 (solid lines) for canalyz- 
ing C2 networks (stars) andCi networks (squares). The curves 
for the different types of networks are plotted in such a way 
that the abscissae are of the order of the number of relevant 
nodes for each type. The data were extracted from those used 
for Fig. 3. The dotted line corresponds to the averaging over 
1000 and not 2000 network realizations. 



a network having typically of the order of not more than 
10 relevant nodes. This explains why for N re i > 10 
the number of found attractors decreases with N re i in 
contrast to the analytical result. On the other hand, the 
analytical result is only valid for N re i 3> L, so that it 
is simply impossible to see the predicted power law with 
computer simulations. A remarkable feature of Fig. 7 
is the qualitative difference between the curves for even 
and odd L. We know from [12] that Gl is smaller for 
odd L than for neighboring even L. 

Fig. 8 shows the total number of attractors found 
starting from 1000 initial states for each network. The 
curves for the different types of networks arc plotted in 
such a way that the abscissae are of the order of the 
number of relevant nodes for each type. 

Just as for the previous figure, with ~ 2 10 random ini- 
tial states the relevant nodes of a network with N re i ~ 10 
assume a large proportion of their possible values, and 
the average number of attractors found is a good esti- 
mate of the real ensemble average (if we average over a 
sufficiently large number of network samples, compare 
the dashed and solid line curves in Fig. 8). For much 
larger networks, it is unlikely that we get the same at- 
tractor twice using only 1000 initial states. Therefore the 
average number of found attractors trivially approaches 
1000, yielding no information about the network dynam- 
ics. 

We compared our results with those of [8], where nu- 
merical simulations were performed to obtain median 
number of attractors. This number is in our system less 



FIG. 8: Average number of attractors obtained with 1000 
initial states and 1000 networks (dashed lines) and 2000 net- 
works (solid lines) for canalyzing C2 networks (stars), RBNs 
(circles) and Ci networks (squares). The curves for the dif- 
ferent types of networks are plotted in such a way that the 
abscissae are of the order of the number of relevant nodes 
for each type. The dotted lines correspond to the curves 
20.6JV?- 5 and 2 o.55^ 1 °- 7 ; where ^ ig the number of nodes 

in Ci networks, i.e. the squared value of the abscissa here. 

than the mean number, and the corresponding curve (not 
shown) lies below our data. For C\ networks, a lower 

bound for the average number of attractors is 2 a6v/ ^, 
see [14]. Our data suggest that the number of attractors 
can well be fitted by 2 aNb with two constants a and b. 
The lower bound as well as the fit for C\ networks are 
plotted in Fig. 8 (dotted curves). 

Originally based on computer simulation of small sys- 
tems, Kauffman had suggested that the mean number of 
attractors increases as V~N. For small N, our data are 
compatible with such a relation. The newer and analyt- 
ical results, see f.i. [12], lead to the exponentially large 
number of attractors, also in agreement with our data. 

The similarity in the form of the three curves in Fig. 
8 confirms our understanding of the dynamics in terms 
of the relevant nodes. In the limit N — > 00 the fraction 
of the relevant nodes with two inputs goes to zero for 
the C2 networks and RBNs and the average number of 
attractors is mainly determined by of the order of In N 
relevant loops. The average number of attractors grows 
exponentially fast with the number of relevant nodes. 



5. SUMMARY 

In this paper, we have shown that canalyzing random 
Boolean networks have a frozen core, the size of which is 
comparable to that of other random Boolean networks. It 
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follows that the attractors of canalyzing networks are de- 
termined by of the order of A 1 / 3 relevant nodes, which 
are connected to of the order of lniV relevant compo- 
nents, most of which are simple loops. 

We have explained how the frozen core arises starting 
from self-freezing loops. Furthermore, we have investi- 
gated the numbers and lengths of attractors. From the 
properties of the relevant components it follows that their 
average numbers increase faster than any power law with 
system size. Although attractors of canalyzing networks 
are on average shorter than those of RBNs, extremely 
long attractors can also arise in canalyzing networks. We 
have shown this by constructing a relevant component 
that has the same attractors as a relevant component of 
an RBN. All the results valid for C2 networks apply also 



to mixed C\ /C2 networks. 

We have also seen that incomplete sampling leads to 
large fluctuations and uncertainties in the data. Addi- 
tionally, very short and very long attractors are difficult 
to find. The first ones constitute an exponentially small 
fraction of the state space, the others appear exponen- 
tially seldom in a network realization. Therefore, com- 
puter simulations needed to be supplemented by analyti- 
cal arguments. At the same time, it is extremely difficult 
to verify numerically some known analytical results. 

We conclude that the original hypothesis that critical 
canalyzing networks have short attractors cannot be up- 
held. Rather, biological systems need to be modeled by 
more specific networks that are not randomly connected. 
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